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Abstract 

The multilayer multiconfiguration time-dependent Hartree (ML-MCTDH) method is discussed 
and a fully general implementation for any number of layers based on the recursive ML-MCTDH 
algorithm given by Manthe [J. Chem. Phys. 128, 164116 (2008)] is presented. The method is 
applied first to a generalized Henon-Heiles (HH) Hamiltonian. For 6D HH the overhead of ML- 
MCTDH makes the method slower than MCTDH, but for 18D HH ML-MCTDH starts to be 
competitive. We report as well 1458D simulations of the HH Hamiltonian using a seven layer 
scheme. The photoabsorption spectrum of pyrazine computed with the 24D Hamiltonian of Raab 
et. al. [J. Chem. Phys. 110, 936 (1999)] provides a realistic molecular test case for the method. 
Quick and small ML-MCTDH calculations needing a fraction of the time and resources of reference 
MCTDH calculations provide already spectra with all the correct features. Accepting slightly larger 
deviations, the calculation can be accelerated to take only 7 minutes. When pushing the method 
towards convergence, results of similar quality than the best available MCTDH benchmark, which 
is based on a wavepacket with 4.6 x 10 7 time-dependent coefficients, are obtained with a much more 
compact wavefunction consisting of only 4.5 x 10 5 coefficients and requiring a shorter computation 
time. 
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I. INTRODUCTION 



The multiconfiguration time- dependent Hartree (MCTDH) method [l-5| has become over 
the last decade the tool of choice to accurately describe the dynamics of complex multidi- 
mensional quantum mechanical systems, providing in many cases reference results that are 
used to benchmark other approaches or determine the behavior of approximate methods. 
MCTDH was initially formulated with molecular quantum dynamics in mind, where the 
degrees of freedom are distinguishable and the potential operator correlates in principle all 
vibrational degrees of freedom of the system. Many successful applications over the years 
deal with the spectroscopy of molecules js _ lQ|, isomerisation and intramolecular vibrational 
energy redistribution (IVR) Ijj], inelastic and reactive scattering calculations 13-17], 



scattering of atoms or molecules at surfaces 
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etc. More recently, the potential of the 



MCTDH method for the treatment of systems of indistinguishable particles, either fermions 
or bosons, has been realized as well, and many applications can be counted nowadays in such 
new directions 2jJ-|27J . Although the equations of motion remain the same in such cases, the 



symmetry properties of the wavefunction and the often two-body nature of the interactions 
has lead to the ^appearance of dedicated programs that specifically and efficiently deal with 
such cases 
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30|. 



In general, the solution of the time-dependent Schrodinger equation in a direct-product 
basis of ID functions (primitive functions) scales exponentially as N* , where N is the number 
of primitive basis functions per degree of freedom and / is the number of degrees of freedom. 
In the MCTDH ansatz one introduces an optimal time-dependent (TD) basis (called SPF 
for single particle function) for each degree of freedom, which can be kept smaller than the 
underlying primitive basis, leading to a better scaling of the number of configurations w . 
In this way, MCTDH can deal with larger systems than the standard method, although the 
exponential scaling with / remains, but to a lower base. The base can be further reduced 
by grouping degrees of freedom together in what are called combined modes or logical 
coordinates. One obtains in this way a smaller number of effective degrees of freedom p, 
leading to a smaller number of configurations, but the TD basis functions that have to be 
now propagated are multidimensional. For a given problem, it is possible to find mode- 
combination schemes that provide an optimum balance between the effort of propagating 
the expansion coefficients of each configuration (the A- vector in usual MCTDH terminology) 
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and the effort of propagating the multidimensional SPFs . With increasing dimensionality, 
however, the multidimensional SPFs, the A-vector, or both, become harder and harder to 
propagate. As a rule of thumb and very generally, the number of degrees of freedom that can 



be treated for a correlated vibrational problem nowadays is about 20, but for systems we 



suited for MCTDH more than 60 degrees of freedom (DOF) have been accounted for |3ll.l32|. 

Multilayer (ML) MCTDH represents a powerful extension of the usual MCTDH approach, 
in which a multiconfigurational ansatz is used for the multidimensional SPFs themselves. 
This results in an extra layer of TD coefficients that have to be propagated and hence 
the name. One way to think of the ML wavefunction ansatz (given below) is to see the 
standard propagation method and MCTDH as particular cases of a ML wavefunction. In 
the standard propagation method a single layer of expansion coefficients with respect to 
some time-independent (TI) basis is present. MCTDH contains a first layer of expansion 
coefficients with respect to the SPF basis, and a second layer of TD expansion coefficients 
that parameterize the time evolution of the SPFs. Expanding multidimensional SPFs using a 
MCTDH-like ansatz one obtains a three-layer scheme. Schemes with more than three layers 
are of course possible. In passing we note that the ML-MCTDH ansatz can be connected 



33] and to density matrix 



to tensor contraction techniques in the mathematical literature 
renormalization group (DMRG) theory 34 1. 

ML-MCTDH was first formulated by Wang and Thoss, who also provided an implementa- 
tion of the method for three layers (in the present paper MCTDH is regarded as a two-layer 
scheme) and showed its applicability on the spin-boson model and on an electron-transfer 



model 35]. A formally identical formulation, but without an implementation, was given 
independently at the same time by Meyer and Worth who termed it cascading MCTDH. 
Over the last years Wang and Thoss have further developed ML-MCTDH, expanding the 
number of layers that their code can handle and treating condensed phase model Hamil- 
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The same 



tonians of large dimensionality between 100 and a few thousand DOF 
authors also described how ML-MCTDH can be used to deal with identical particles, which 
is accomplished defining the problem in Fock space and directly working in the occupation- 



number representation 



40j. (Note that ML-MCTDH cannot be applied to identical particles 



in the usual first-quantized formulation because the grouping of DOF into combined modes 
destroys the necessary exchange symmetry properties of the wavefunction.) Recently, Man- 
the provided a recursive formulation of the ML-MCTDH equations of motion (EOM) for 
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any number of layers and described the recursive algorithm that has to be used to compute 



all quantities entering the EOM 
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42|. 



In this contribution, we present the implementation of the EOM and recursive algorithm 
proposed by Manthe, and discuss its applicability and performance on a Henon-Heiles (HH) 
model system and on a realistic non-adiabatic molecular case, full dimensional 24D pyrazine. 
The resulting implementation is fully general and it can handle standard method (one layer), 
normal MCTDH (two layers) and any desired multilayering scheme. In this paper, we 
report on simulations with three layers and up to a seven-layer case for the 1458D HH 
Hamiltonian. We also simulate pyrazine using five and six layer schemes. The ML-MCTDH 
concept, EOM and recursive algorithm are discussed first, emphasizing some aspects of our 
implementation such as the treatment of the separable parts of the Hamiltonian. We have 
tried to keep such discussion self-contained and to reflect the way our code operates, but 
for an in-deep discussion Ref. 4l| should be read. The performance and correctness of the 
recursive implementation are tested on the HH Hamiltonian. For this system, 6D calculations 
with varying coupling strength are carried on first, which are directly compared to MCTDH 
results. For such a low dimensionality, the algorithmic complexity of ML-MCTDH does not 
yet pay off, and MCTDH is more efficient than ML-MCTDH. For 18D HH we find that ML- 
MCTDH starts to be competitive. For this case we analyze the convergence properties of the 
method with respect to the number of basis functions at each layer. Simulations of 1458D HH 
are also reported, showing the power of this approach. The ML-MCTDH implementation is 
then applied to the photoabsorption spectrum of pyrazine using the 24D vibronic-coupling 
Hamiltonian of Raab and coworkers 6| . This is a realistic molecular Hamiltonian involving 
the presence of a conical intersection between the electronic state S2 initially populated by 
photoabsorption and the Sj electronic state to which the system decays within some ten fs. 
The Hamiltonian of Ref. |6| has been used to test several quantum dynamical approaches. 
In the case of ML-MCTDH, we find that the method performs very efficiently and provides 
already well converged spectra in about a fifth of the CPU time used by the smallest of the 
MCTDH reference results using extremely few computational resources. Tolerating a larger 
error one can reduce the computation time to 7 minutes and still get a spectrum with all 
basic features in place. This is amazingly fast, considering that one is dealing with a fully 
correlated 24-dimensional quantum dynamics calculation. No other quantum dynamical 
method applied up to date to the pyrazine Hamiltonian offers the extreme quality/cost 
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relation provided by ML-MCTDH. When pushing ML-MCTDH towards convergence with 
the pyrazine Hamiltonian, we find that it still offers results of at least the same quality (if 
not better) than the best benchmark results available from MCTDH calculations with about 
4 x 10 7 configurations, at a lower computational cost. 

The paper is organized as follows: Section [III discusses the ML-MCTDH equations and 
implementation. Section IIII Al discusses the results on the HH Hamiltonian and Section 
IIIIBI discusses the application of ML-MCTDH to the photoabsorption spectrum of pyrazine. 
Section [TV] summarizes the results and provides some future perspectives for the method. 



II. MULTILAYER MCTDH 

A. Ansatz and general concept 

The wavefunction of a system of distinguishable degrees of freedom, e.g. representing 
molecular vibrations, can be conveniently represented as an expansion in terms of direct 
products of orthonormal basis functions, one for each degree of freedom 

tffa, . . . , q f , t) = E ■ ■ ■ E 4 (0 • Xfffoi) • • • • • Xf f \q f )- (1) 
h if 



Following the notation introduced by Manthe 



411 ] . the Aj^j denote the TD coefficients of 



the first (and by now the only) layer of TD expansion terms. 

In the MCTDH approach, the wavefunction is expanded in terms of direct products of 
orthonormal time- dependent single particle functions (SPF) 

tffa, . . . ,q f ,t) = E " i/W • V^ ;1) (<?!,*) • • • • • ^ J \q f ,t), (2) 

h if 

which are themselves expanded in terms of the underlying primitive basis 

^(?«^) = E A S(*)-4 K) (^)- ( 3 ) 

Therefore, MCTDH can be seen as a 2-layer scheme with TD coefficients A 1 ^ ^ (t) at the 
top layer, and sets of second layer TD coefficients A^At) for each degree of freedom. We 
usually refer to the 1-layer scheme as the standard method, to the 2-layer scheme simply as 
MCTDH, and to deeper layering schemes as ML-MCTDH. Note the important detail that 



all SPFs of a certain degree of freedom are expanded in terms of the same underlying basis, 
i.e. the Xj (in) functions in Eq. ([3]) have no m index. 

The computational gain of MCTDH with respect to the standard method arises from the 
expansion orders n K being in general smaller than the size of the underlying primitive basis 
N K , which leads to a smaller number of TD coefficients to be propagated. However, the 
total number of TD coefficients is still given by n«=i n <« an d therefore the computational 
effort raises exponentially with the number of degrees of freedom. Hence, MCTDH does 
not eliminate the exponential scaling but reduces the base to which the scaling occurs (for 
a detailed analysis of the computational scaling of MCTDH see for example Ref. 3|). The 
base can be further reduced by combining the physical coordinates q\,...,qf into logical 
coordinates (also referred to as combined modes) Q\, . . . ,Qp, such that each logical coor- 
dinate comprises one or several of the physical coordinates, as Q\ = {Qi K , ■ ■ ■ ,Qd K }- Here 
a more general notation to designate a combined coordinate has been introduced, that will 
be useful when discussing deeper layering schemes. Mode combination has been extensively 
used in many applications of the MCTDH method. Similarly as in Eq. ([2]), the MCTDH 
wavefunction now reads 

*(Ql, ■ ■ ■ ,Ql,t) = £• ■ ,.(/) • ^\QU) • . . . ■ ^ p \Ql,t), (4) 

and the TD basis functions ipj^ K \Q\, t) are now multidimensional. Introducing mode- 
combination, the computational effort is switched from the propagation of a large vector 
of A] „• (t) coefficients and one-dimensional SPFs, to a shorter vector of coefficients but 
multi-dimensional, harder to propagate SPFs. It is often the case, that some experience 
and knowledge of the system under consideration is necessary to find an efficient mode- 
combination scheme for a given problem. The mode-combined SPFs are given by 

Ni K N dK 
h K 3d K 

where the last equation resembles the ansatz in Eq. (JTJ) in that multidimensional SPFs are 
represented as a multiconfigurational expansion in terms of underlying time-independent 
basis functions. In the same way a TD basis was introduced above in going from the 
standard ansatz to MCTDH, effectively adding a second layer of TD expansion coefficients, 
an MCTDH expansion can be used to represent the ipl£{Q\,t) SPFs, effectively adding a 
third layer of expansion coefficients. This results in a 3-layer ML-MCTDH ansatz. 
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The ML-MCTDH layering scheme can be very flexible. In a high dimensional system one 
combines groups of degrees of freedom into high dimensional SPFs until the size of the vector 
of coefficients in Eq. (j3J) becomes manageable. However, the combined SPFs are too large to 
be efficiently propagated. Then, one breaks the combined modes into even smaller groups 
of logical coordinates introducing a new layer of coefficients, whose size is manageable. This 
procedure can be repeated over and over until the (possibly combined) primitive degrees of 
freedom are reached. 

Owing to the flexibility of the layering scheme, and the fact that ML-MCTDH wave- 
functions can be many layers deep, it is convenient to introduce a diagrammatic notation 



to represent these objects 



411 ] . In this notation, wavefunctions are represented by trees, 



i.e. connected graphs without loops. Each node in the tree represents a set of vectors of 
coefficients ^^^".'.'.j" 1 , where I denotes the layer depth, Ki, . . . , ki_i denote the indices of 
the logical degrees of freedom starting from the top node and down to a particular node 
(path from the top down to a given node), m indicates the different vectors within this node, 
and finally ji, ■ ■ ■ ,j PK are the tensor indices of each particular array of coefficients within 
the node. 

Some ML diagrams are presented later when discussing the various applications (see 
Figs. [HUH])- The lines connecting one node with its descendant nodes in such diagrams 
represent the tensor indices jx, ■ ■ ■ ,jp K1 K , one line per index, and the numbers at the side 
of each line refer to the maximum possible value of the corresponding index. Each node is 
uniquely described by the values of its label (/; k±, . . . , k;_i). 

B. ML-MCTDH equations of motion and recursive implementation 

The ML-MCTDH equations of motion (EOM) have been derived and discussed by Wang 



and Thoss [35| and Manthe |4l|. In Ref. |4l| a fully general derivation of the ML-MCTDH 
EOM for arbitrary layering schemes is provided, together with an algorithm for the recursive 
evaluation of all intermediate quantities entering the EOM. We reproduce here the EOM 
and hint the key elements of the recursive algorithm for the sake of completeness. Later we 
discuss some specific aspects of our implementation. 

The ML-MCTDH EOM have a very similar structure to the usual MCTDH equations, 
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and for the top layer coefficients they are identical to the MCTDH ones 

~dt 



= J>}|£|$*)^, (6) 



where the top layer configurations $j = t Pj 1 1 ' 1 \Q\,t) ■ . . . ■ ip^ p ' p \Qp,t), are defined as direct 
products of SPFs and the multi-index J = ji, ■ ■ ■ ,j p has been implicitly introduced. The 
EOM for the propagation of the SPFs are formally the same for all layers 

= (1 - P$ 5> 2A ' W " (HVj'X 1 , (7) 

j,m 

where = '%2j\ ( Pj' Kl ){ < Pj' Kl \ * s the projector onto the space spanned by the (p* ,Kl SPFs, 
p z,Kl is a density matrix and (H) z,Kl is a matrix of mean- field operators acting on the ipj' Kl 
functions. The symbol z is a shorthand for the indices specifying a node, z — l; K\, . . . , Ki-i, 
and for further reference we introduce z — 1 = I — 1; K\, . . . , Ki_ 2 , and similar for z + 1. In 
its form above, the EOM for the SPFs look identical to the usual EOM for the SPFs in the 
usual MCTDH 3]. Only the computation of the density matrices and mean-fields entering 
the EOM is now more involved than in a single-layer MCTDH scheme. 

Before describing how in practice the recursive calculation of mean-fields and density 
matrices is performed, we first restrict ourselves to the case in which the Hamiltonian is 
given by sums of products of operators acting on the primitive (but possibly combined) 
degrees of freedom Q^ rim 

s s f 

H = °rHr = °r W ' ( 8 ) 

r=l r=l ft=l 

This form of the operator is necessary to avoid the computation of high dimensional integrals, 
and has been thoughtfully discussed in the MCTDH literature. Terms which usually do not 
follow the product form are, e.g. general potential energy surfaces. There exist however 



algorithms to bring general potentials to product form [43M46l| . Furthermore, MCTDH and 
ML-MCTDH formulations based on time-depende nt g rids that by-pass this restriction (so- 
called CDVR approach) have also been provided 4l|, |47| , but here we always will assume 
Hamiltonians in product form, Eq. (JS]). 

In order to discuss the practical computation of the mean-field elements (H) '^''"' Kl ap- 
pearing in Eq. ([7]), we center our attention on a particular node z — (I; Ki, . . . , This 
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node relates to the logical coordinate Q K ilf' 1 ''"' Kl 2 , which can be further decomposed into a 
combination of logical coordinates 

= {<&...,<%„}, (9) 

and holds vectors of coefficients denoted by A z m .-^ - p . These coefficients are expansion 
coefficients of SPFs in layer / — 1 expanded in terms of SPFs in layer /. This is better seen 
by extending Eqs. (13|5|) to the fully general multilayer case 

^>^{QtL) = e--e^i...^ n ( ioa ) 

= £^-*M£i)> (iob) 
j 

where the configurations $5 are the multilayer generalization of the configurations intro- 
duced in Eq. (jBJ). The operator summands in Eq. (JSJ) can be written as 

H r = n^ u -' Kl ■ hi* 1 '-'" 1 (lla) 
P K 1+1 

= u^ 1 '-'* 1 ■ jj ( llb ) 

«i+l=l 

where h l j Ku '" ,Kl acts on the logical coordinate Q 1 k* 1 '"'' Ki ~ 1 , while % l j K i'-' K i acts on the corre- 
sponding remaining space. As seen from Eq. (Ill b ). h^ K1 ''"' Kl can as well be written as a 
direct product of operators hr^ 1 '* 1 ''"'* 1 * 1 , which act on logical coordinates Q l + l '^-^ K i _ i n 
fact, both operators i-l l J K ^>-> K i and h l j KU '"' Kl can be broken down to a simple product of oper- 
ators which act on one primitive (but possibly combined) coordinate only (see Eq. (JSJ)). The 
matrix elements of the mean-field operators at the right-hand-side of Eq. [7] for one particular 
summand H r are given by 

(Hr)ffij=(^\(Hr)'£\*?), (12) 

which is readily seen if the SPFs in Eq. ([7]) are written in their explicit form given in 
Eq. (110b ). These are the most cumbersome quantities to compute within the ML-MCTDH 
scheme. By substituting Eqs. (ITTk.b) into Eq. (fl2|) one arrives at 

(H r )t)iJ Kl = ^r- R i£- Kl ■ [/''r _ /; (13a) 

Pk I+1 

= (ii^ i '-' Ki ) l j^'---' Ki ■ Y[ \h l r +1 ' Ki '---^+A (13b) 

« I+ i=i ^+1.^+1 



with 



I J 



($f>^l$* +1 ) 



(14a) 
(14b) 



1 *1+1,Jk i+1 

(The mean-field tensor ("H^ ,K ')^ ! will be discussed below, Eq. (1TB]) ). It remains therefore to 
be seen, how the quantities appearing in the right-hand-side of Eq. (113b ) can be recursively 
computed. 

The expression for the /i-matrices is found by starting from Eq. (114b) and explicitly 
writing the SPFs in layer / + 1 in terms of those in layer I + 2 



{<P„ 

\ra K , 



Z+1,K/ + 1 I 



i+1 



n ^ 2 

n EE- 4 '*' 



^+2|, Z+1 ' K '+l\ 



(15) 



^+2 



; i+2 



.4 



--+2 



i+2 • 



^(+2=1 J K; + 2 i,j 

Here multi-indices J K and have been introduced to make the notation more compact. 
The former contains all indices except for the K-th one, while the latter has index k set to 
m. Eq. (j!5p tells us that the computation of the matrix elements of hr +1 ' Kl+1 requires the 
previous knowledge of the matrix elements of hr +2 ' K ' +2 . Therefore, one starts at the lowest 
layer with the matrix elements of H with respect to the underlying time-independent basis 



h 



<xl K Vr m;K ixr> 



(16) 



where Z max denotes the maximal layer, i.e. Z max = 2 for normal MCTDH, and proceeds 
bottom up through the tree structure building the matrix representations of the ^ +1,reir "' Ki + 1 
operators, which will subsequently be used to build the matrix elements of h l j KX >-< K i operators, 
and so on. 

The mean-field matrix elements are computed in practice in a similar way to Eq. (fT5"j) . 
For the top layer, the mean-field elements are given as 



h 



h3 



A 1 



(17) 



where the composite indices J K,U and J^ n are defined similarly as J K and J^. Eq. (117)1 
corresponds to the normal way of calculating the mean-fields in standard MCTDH, and is 
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an alternate way of writing, e.g. Eq. (67) in Ref. [3|. For a given deeper layer, the calculation 



.a. 



involves the mean-fields of the above layer and is given by 



4l| 



Z,Kl 

m,n 



a,b 



A* 

/ -> a;J, 



h 



2— 1,1/ 



.4 



1,3 



(18) 



where a and b run over the (I — l)th layer SPFs, which are contained in the Ith layer 
SPF m and n, respectively. Therefore, the mean-fields can be computed if the mean-fields 
of the layers above, as well as the matrix elements of the h operators, are known. The 
density matrices are computed in a very similar way to Eq. (fT8"|) . but application of the h 
matrices to A vectors on the right side of Eq. (IPS]) disappears. Eqs. (I17H18P are obtained 
after introducing the so called single-hole functions, which are defined as the variation of 
the total wavefunction with respect to a particular SPF. We leave this more specialized step 
out of this discussion and refer the interested reader, e.g. to Refs. j^| 



or 



41|. 



The recursive algorithm proceeds then as follows: 1) Starting from the bottom of the tree 
and proceeding up, the h- matrix elements are evaluated after Eq. (I15p . 2) Starting from the 
top of the tree and proceeding down the mean-fields and density matrices are computed as 
in Eq. [18j 3) The right-hand-side of the EOMs, Eqs. f[6"fT|) . is evaluated. 



C. Separable part of the Hamiltonian across layers 

The efficiency of the ML-MCTDH calculations is increased if one takes advantage of the 
terms fif"™'^ in Eq. flH]) that are unit operators. The kinetic energy part of the Hamiltonian 
contains usually many unit terms. Also common models such as Henon-Heiles, vibronic 
Coupling, system-bath, etc. have as well many unit terms in the potential energy part of 
the operator. 

The obvious way to obtain computational savings is to flag all unit terms as such, so 
that they do not have to be stored as unit matrices and they do not have to be explicitly 
multiplied. This is done from the bottom of the tree and proceeding recursively upwards. 
If all hr +1 ' Kl+1 = 1, then hp Kl is flagged as 1. Our implementation keeps track of all unit 
fo-terms at all layers, which saves resources in the recursive construction of the /i-matrices 
in Eq. ([To]) , in the mean- fields construction in Eq. ( 1181) and in the equations of motion. 
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Once the unit operators at all layers have been flagged as such, it is possible to identify 
non-correlated terms for a particular node, i.e. terms that decompose as 

H r = l z ' Kl ■ h z r ' K \ (19) 

such that the non-unit parts of H r operate only on the combined coordinate Q Z K[ . The step 
of identifying separable terms is conveniently done from the top of the tree and proceeding 
downwards. To illustrate this let us refer to a simple 4D example and the operator summand 

H r = h^' 1 ■ if im ' 2 ■ h p / im ' 3 ■ /*P rim;4 . (20) 

The considered tree structure has four primitive nodes (3; 1,1), (3; 1, 2), (3; 2, 1) and (3; 2, 2) 
with /i-terms hf' 1,1 , l^' 1 ' 2 , hy 21 and hp 2 ' 2 , which are the matrix representation of the primitive 
product terms in Eq. (|2"U1) . The two nodes at the intermediate layer, namely (2; 1) and (2; 2), 
have therefore fo-terms h 2 ' 1 ' 1 and l 2 ' 1,2 for node (2; 1) and /i-terms h 2 ' 2 ' 1 and h 2 ' 2,2 for node 
(2; 2). The top node (1; ) has /i-terms h]} 1 and hi' 2 . Only in the case that both h 2 ' 1 ' 1 and 
l 2 : 1 ' 2 would have been unit operators, would hi' 1 have been marked as a unit operator as 
well. Now we turn to the identification of non-correlated terms starting from the top of 
the tree. H r is marked as correlated in node (1; ) because both hp 1 and hi' 2 are different 
from the unit operator. At first sight, the H r term could be flagged as non-correlated in 
node (2; 1) because only one of the /z-terms is different from unity. However, the parent 
node had this term marked as correlated, and therefore H r is marked as correlated in node 
(2; 1). The reason why this happens is that H r acts also on node (2; 2), which is in another 
branch of the tree. It is then clear that a convenient way to spot situations of this kind 
is to proceed recursively from the top of the tree, as illustrated. Ultimately, the H r term 
is marked correlated in node (2; 1) because the mean-field matrix elements (^v' 1 ) 2 ^ 1 and 
("H 2 ' 1 ' 2 ) 2 ^ 2 correlate the motion of SPFs pf 1 ' 1 and y? 2,1,2 with the dynamics proceeding in 
node (2; 2) and below. The coupling to the dynamics in the other branch arises from the 
relation of the mean-field elements just mentioned to the mean-fields {7il' ,1 )l^ n in the top 
node, as shown in Eq. (1TB]) . 

Once the non-correlated terms have been identified for each node, one can exploit, in the 
same way as in standard MCTDH, that for them the mean-field matrix elements become 
identical to the elements of the density matrix 

(¥wr--,«iym,.,*i = (^«.-.«i) mn . (21) 
12 



Hence, for a particular node with H r terms 1 to s non-correlated, and terms s + 1 to t 
correlated, Eq. ([7]) can be rewritten as 



dt 



\g=l j,m r=s+l / 

which is the EOM for the SPFs on which our implementation is based. The EOM for the 
coefficients follows directly from this equation and Eq. (ITU]) 

There is no contribution from d<&j/dt because, due to the projector, the SPFs are orthogonal 
to their time derivatives. Hence one may write compactly 

dA z 

-g 1 = E ~ E • , J,/- , ;:a ) E KmKL K,;L , (24) 

K j m,L 

where a detailed expression for the matrix M^ m . JL follows from Eqs. (|2"2l [231 HU]) . 

As a final comment, it is worth mentioning that in large systems with groups of coor- 
dinates strongly correlated among them but weakly correlated to other groups, i.e. with 
some sort of locality, the identification of the correlated terms node by node may become 
very important. In such cases there will be many correlated terms in nodes at lower layers, 
but less of them as one proceeds up to the top of the tree. Such an example is given by be 
the HH model to be introduced below, where each harmonic oscillator is coupled only to 
its next neighbors. For similar reasons, when dealing with system-bath Hamiltonians with 
ML-MCTDH, it may be advantageous in most cases to separate system and bath already 
at the top layer. 



III. RESULTS AND DISCUSSION 



The ML-MCTDH algorithm and implementation described above are applied to two 
paradigmatic systems. On the one hand we simulate the HH model for various dimension- 
alities, and on the other hand we test our implementation on the pyrazine system using the 
second-order vibronic-coupling Hamiltonian of Ref . jsj . 
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A. Henon-Heiles 



The HH Hamiltonian 

k = f E ("^2 + + A E " 3&i. J (25) 

written here in dimensionless units, offers a convenient playground to investigate the con 
vergence properties and performance of ML-MCTDH. The HH model was used in Ref. [48 
to benchmark the MCTDH method and the same Hamiltonian with similar parameters is 
used here for the bechmarking of ML-MCTDH. For all simulations we set u — 1. The degree 
of anharmonicity and coupling are controlled by the single parameter A, for which different 
values are chosen. The initial amount of energy in the system is easily controlled by the 
position at which the initial wavepacket is centered for each degree of freedom. The initial 
wavepacket is always taken as a Hartree product of Gaussian functions of width correspond- 
ing to the ground vibrational state of the harmonic part of the Hamiltonian. Initially, all 
momenta are set to zero, (p K ) = 0, and similarly all positions, except for a few coordinates 
which are displaced by two length units with respect to q K = 0. The HH Hamiltonian in 
Eq. f )25|) has dissociative channels due to the cubic terms, which are accessible in the energy 
range of the reported simulations. To avoid reflections from the grid edges, the Hamiltonian 
is augmented with a complex absorbing potential (CAP) as H = T + V — iW. Further details 



on the CAP used are found in Ref. 



48 



]. A sine-DVR (for discrete variable representation) 
primitive basis with points between —9 and 7 length units and A = 24 grid points is used for 
all coordinates throughout the HH simulations. All propagations on HH models are carried 
on up to 30 time units, which roughly corresponds to five oscillation periods T = 2ir/u in 
the harmonic limit. 



1. 6D simulations 

FIGURE IB AROUND HERE 

We first concentrate on a 6D HH model. Calculations with A = Ao and A = 2Ao are 
performed, for Aq = 0.111803. A value of Aq = 0.111803 allows to compare with other works 



that use the same Hamiltonian 



48j. For the reported 6D calculations, the initial wavepacket 
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is centered at q K = 2 for coordinates q 2 , q± and q§. Standard MCTDH (2- layer) and ML- 
MCTDH calculations using a 3-layer scheme are compared, and the trees representing the 
MCTDH and ML-MCTDH wavefunctions are displayed in Figs. [T^, and [lb, respectively. 
In the MCTDH case every two primitive DOF are combined, and N 2 equals the number of 
DVR grid points (primitive basis functions) for each DOF. In the ML-MCTDH wavefunction 
an extra layer of TD coefficients is introduced, effectively representing the 2D SPFs of 
the combined modes in the original MCTDH wavefunction by a new multiconfigurational 
expansion. It is worth noting that in the case of having N 2 = N 3 for the ML-MCTDH 
wavefunction, the lowest layer is treated exactly and both MCTDH with mode combination 
and ML-MCTDH provide numerically identical results (for the same number of SPFs JVi at 
the top layer). 

The different propagations are compared by inspecting the autocorrelation function ob- 
tained as a(t) = (ty*(t/2)\ty(t/2)). (This trick j^] yields an autocorrelation function which 
is twice as long as the propagation). The autocorrelation at long times is difficult to con- 
verge, since the (ML-)MCTDH wavefunction accumulates error as time increases, providing 
an adequate quantity to compare the different propagations. Conversely, the spectra cor- 
responding to the autocorrelation functions have the error averaged over the whole energy 
range, and therefore are not such a good quantity to look at for a stringent comparison. 

FIGURE M AROUND HERE 

FIGURE SI AROUND HERE 

Results for the low (A = Ao) and high (A = 2Ao) coupling regimes are shown in Figs. [2] 
and El respectively, where the absolute value of a(t) is displayed. In both Figs. [2] and [31 the 
numbers in parenthesis correspond to (N%, N 2 ) in the trees in Fig. [TJ For MCTDH N 2 equals 
the number of grid points and is given only for completeness. In the case of low coupling 
(Fig. |2J), all calculations yield almost indistinguishable autocorrelation functions. Fig. [2b 
displays the difference between autocorrelation functions of the ML-MCTDH calculations 
and the best converged (50,24) MCTDH calculation. The (50,20) ML-MCTDH calculation 
is very close to the reference MCTDH result. The calculations that diverge earlier with 
respect to the reference result are the (30,10) and (50,10) ones due a poorer convergence 
of the lower layer, while the ML-MCTDH (30,20) calculation starts to deviate from the 
reference result after the system has completed five to six oscillation cycles. 
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After doubling the strength of the coupling parameter, the situation becomes more com- 
plex. In Fig. [3^ one sees that after two oscillation periods the autocorrelation functions of 
the different propagations start to diverge on an appreciable scale. Only the best MCTDH 
calculation (50,24), and the best ML-MCTDH (50,20) result in very close autocorrelations 
up to 60 time units. Fig. [3b shows the difference between the ML-MCTDH results and the 
(50,24) MCTDH calculation. The two worst calculations are clearly the ones with N 2 = 10. 
In this case, the correlation at the lowest layer, i.e. between DOF (qi, q 2 ), (#3, #4) and (g 4 , g 5 ) 
is not well represented. It is important to note that the (50,10) calculation does not offer an 
improvement over the (30,10) one. The missing correlation within the logical coordinates 
in the lowest layer cannot be regained by having more SPFs in higher layers. Similarly as 
in standard MCTDH calculations, the eigenvalues of the density matrices at each layer, the 
natural populations (NP), provide an indication of the degree of convergence of a calcula- 
tion. As a rule of thumb, the smallest NP should be of the order of 10~ 5 to provide good 
results in many cases, although this can vary depending on what quantity one is trying to 
compute. For the (50,10) and (50,20) ML-MCTDH calculations at t = 60, the smallest NP 
of the three logical coordinates at the top layer are (1.4 • 10~ 5 / 1.3 • 10~ 4 / 9.3 • 10~ 5 ) and 
(3.4 • 10~ 5 / 1.8 ■ 10~ 4 / 1.3 ■ 10 -4 ), respectively. Therefore, the convergence at the top layer 
is slightly worse for the (50,20) case than for the (50,10) case, even though by comparison 
to the standard MCTDH calculation it is clear that the (50,20) calculation provides better 
converged results than the (50,10). Hence, by allowing more SPFs at the second layer, the 
evolution of the top layer SPFs becomes more complex since they have more variational 
freedom. 

As a remark to the 6D simulations, in this case standard MCTDH calculations are more 
efficient than ML-MCTDH ones. Only in the case of having a larger number of DVR 
points per degree of freedom resulting in very large 2D combined coordinates, would the 
introduction of a further layer (compare Figs.QJi and [lb) constitute an advantage. Efficiency 
issues are discussed later in relation to the larger HH simulations and to pyrazine. 

2. 18D simulations 

FIGURE H AROUND HERE 
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Next we discuss simulations on a 18D HH model with A = Ao- The initial wavepacket 
is centered at q K = 2 for coordinates g 4 , g 8 , q 12 and gi 6 , resulting in an initial energy of 
15.807. As in the 6D case, wavepackets are propagated up to 30 time units. For this model 
we run standard MCTDH calculations in which the 18 physical coordinates are grouped 
into 6 combined modes (logical coordinates) with three primitive DOFs each. The tree 
representation of the MCTDH wavefunction is given in Fig. UK There, N\ is the number 
of SPFs per combined mode and N 2 = 24 is the number of DVR functions per DOF. Two 
standard MCTDH calculations were performed, using Ni = 10 and Ni = 14 SPFs per 
combined mode, resulting in 10 6 and 7.5 ■ 10 6 configurations, respectively. ML-MCTDH 
calculations based on the tree structure in Fig. Hb and using the same Hamiltonian and 
initial conditions were performed as well. In the ML calculations various basis sizes were 
employed, always keeping Ni equal to N 2 . The autocorrelation functions resulting from the 
MCTDH and ML-MCTDH calculations are presented in Fig. [5J where the numbers in the 
legend correspond to JVi for MCTDH, and to Ni and N 2 for ML-MCTDH. At the scale 
shown in Fig. [5^ all simulations yield very similar results. Fig. [5b shows the last revival of 
\a(t) \ in more detail. The worst result appears to be the standard MCTDH calculation with 
Ni = 10. The various ML-MCTDH calculations approach the (Ni,N 2 ) = 20 result as the 
size of the basis increases, while the standard MCTDH with Ni = 14 seems not to be yet 
be fully converged with respect to the trend of the ML-MCTDH simulations. 

FIGURE E AROUND HERE 

The set of simulations described above were repeated doubling coupling parameter 
A = 2Ao- For this coupling strength the revivals of the autocorrelation function decay af- 
ter a few oscillations. The different simulations yield still similar autocorrelation functions 
(Fig. [6^), but differences are now larger than in the small coupling parameter case, which 
can be seen by inspecting Fig. Eb. There, the autocorrelation function between 10 and 30 
time units is shown in detail, and after this time the autocorrelation has almost vanished. 
By looking at the autocorrelation functions of the various MCTDH and ML-MCTDH sim- 
ulations one sees that full convergence with respect to the number of SPFs has not been 
completely reached, even at short propagation times. Full convergence of the 18D HH model 
with a coupling parameter of A = 2Ao is a very demanding task. The N\ = 14 MCTDH 
simulation was run on 8 processors using shared-memory parallelization, for which a speed 

17 



factor of about 3 can be expected |62j. Such calculation required 71 hours of wall clock time 
on the 8 processors and would have taken about 200 hours on a single processor. About 
5 GB of main memory were used for the wavepacket propagation. In contrast to this, the 
largest ML-MCTDH calculation with (A^A^) = 20 required 136 hours of wall-clock time 
running on a single processor and the wavepacket propagation used less than 500 MB of 
main memory, while the (A 1; A 2 ) = 12 ML-MCTDH calculation needed only 7 hours on 
a single processor. All calculations reported in this paper are performed on an quad-core 
Opteron, processor type 2384, 2.7 GHz. 

FIGURE M AROUND HERE 

If we now look at the spectra of the propagated wavepackets obtained by Fourier Trans- 
form (FT) of the corresponding autocorrelation functions we see that they are very similar. 
The spectra of two MCTDH calculations with N± — 6 and Ni = 14 are shown in Fig. [71 
The N\ = 10 case (not shown) yields a spectrum very similar to the larger MCTDH calcu- 
lation. The spectra of two ML-MCTDH calculations are also presented in Fig. [7J namely 
the smallest and largest calculations of the ML series. The spectrum of the MCTDH cal- 
culation with Ni = 6 yields somewhat different intensities in the higher energy range than 
the other simulations. When closely inspecting the lower energy range (Fig. [7b) one finds a 
non-negligible energy shift of the spectrum of the small MCTDH calculation to higher en- 
ergies. Interestingly, the Ni = 6 standard MCTDH propagation required 17 hours running 
on 8 parallel processors in the same conditions as discussed above. The number of SPFs for 
this simulation was chosen such that the cost would be comparable to that of the smaller 
(N U N 2 ) = 12 ML-MCTDH calculation. The smaller ML-MCTDH calculation required 7 
hours of wall-clock time on a single CPU running on the same hardware, and its spectrum 
is already very similar to the spectra of the larger MCTDH and ML-MCTDH calculations. 

FIGURE m AROUND HERE 

The size of the 18D HH model is such that standard MCTDH calculations can still be 
conducted to a good accuracy. Contrary to the 6D case however, for this dimensionality it 
already pays off to introduce a further layer in the calculation, and 3-layer ML calculations 
offer results of similar quality than the standard MCTDH ones with a noticeably smaller 
cost. 
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3. 1458D simulations 

Finally, we report ML-MCTDH simulations on a HH model with 1458 DOF and coupling 
constant A = Ao- The system is described in this case by a 7-layer wavefunction. At the top 
level the coordinates are divided in three groups of 486 coordinates, which are subsequently 
divided again in three groups, and so on. This is repeated until groups of two primitive 
DOF are reached in the seventh layer, which are then kept combined. For this deeply 
layered wavefunction two limiting cases are numerically investigated. In example (1), DOF 
§486 and #487 are displaced for t = such that (q K ) = 2. Since the two degrees of freedom 
are adjacent, there is a coupling term in the Hamiltonian between them. However, and 
^487 belong to different logical coordinates at all layers of the tree. They belong to logical 
coordinates Q\ and Q\ at the top layer and to coordinates Q3' 1 ' 3 ' 3 ' 3 ' 3 and Q{ 2, ' ' ' at the 
bottom layer, respectively. In example (2), DOF q 72 g and g 730 are initially set at (q K ) = 2. 
Again, these two DOF are coupled in the Hamiltonian, but in this case they belong to the 
same logical coordinate at all levels, starting from logical coordinate Q\ at the top layer and 
down to the same combined mode Q2' 2 ' 2 ' 2 ' 2 ' 2 a t the bottom layer. In both cases, decoupled 
reference calculations are also conducted in which the coupling term between DOF q^Q 
and ^487 for case (1) and between DOF g 72 g and §730 for case (2) are eliminated from the 
Hamiltonian. The simulated system corresponds to a long chain of coupled oscillators. In 
both cases (1) and (2), the initially displaced degrees of freedom are far from the ends of 
the chain. Therefore the system dynamics in both cases should be very similar, since the 
extremes of the chain will play a role only at much longer times than simulated. 

FIGURE E AROUND HERE 

Fig. Ek, shows the autocorrelation function of the two decoupled reference calculations. 
In both cases the wavefunction consists of 5 SPFs for each logical coordinate at all layers, 
which results in 2 326 520 TD coefficients. Since there is no coupling term in the Hamiltonian 
between the two initially displaced DOF, their distance along the tree structure should be 
irrelevant, and in fact both simulations yield nearly identical results. When the interaction 
is turned on, however, the distance between two coupled DOF along the tree does matter. 
In case (2) the tree distance between DOF g 72 g and q 730 is the shortest possible. Fig. [Hb 
shows the autocorrelation of two propagations for case (2). One of them is based on the same 
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wavefunction as the reference results. A second propagation is based on a wavefunction with 
20 SPFs for the combined mode containing q 72 g and g 730 , Q2' 2 ' 2 ' 2 ' 2 ' 2 , and 15 SPFs for the two 

6"2 2 2 2 2 6'2 2 2 2 2 

neighboring combined modes Q{ ' ' ' ' and Q 3 ' ' ' ' ' . One layer above the wavefunction has 
10 SPFs for logical coordinates Ql' 2,2 ' 2,2 and its two neighbors Q^' 2 ' 2 ' 2 ' 2 and Q3' 2 ' 2 ' 2 ' 2 . Further 
up in the tree, logical coordinates containing ^729 and q^o and its two neighbor coordinates 
have 7 SPFs each. All other logical coordinates at all layers are represented by 4 SPFs. This 
scheme results in 1 853 310 TD coefficients. As seen in Fig. [8b, both wavefunctions yield 
very similar autocorrelations until the end of the simulation. By inspecting the natural 
populations at the final time for both case (2) simulations one notices however that the 
convergence at the lowest layer is rather poor for the homogeneous wavefunction with 5 
SPFs overall, with a population of the order of 10~ 2 for the last natural orbital of coordinate 

6*2 2 2 2 2 

Q2 ' ' ' ' . When inspecting the populations of the propagation with more SPFs for logical 
coordinates close to the displaced mode Q2' 2 ' 2 ' 2 ' 2 ' 2 , substantially improved values of at least 
10~ 4 are obtained for all logical coordinates containing Q2' 2 ' 2 ' 2 ' 2 ' 2 at all levels. 

Example (1), in which the distance between both initially displaced DOF 5436 and 5487 is 
the largest possible within the tree, has a different convergence behavior. Fig. |8t compares 
two autocorrelations for case (1). A first propagation is based in the same homogeneous 
wavefunction already discussed for case (2). A second simulation is based on a wavefunction 
with 15 SPFs for logical coordinates Q^' 1 ' 3 ' 3 ' 3 ' 3 and Q{ 2 ' ' ' 1 ' , which respectively contain 
DOF ^486 and ^437, and with 7 SPFs for logical coordinates containing also g486 and ^437 at 
all layers. All other logical coordinates have 4 SPFs and total number of TD coefficients 
is now 1810 940. Now, both autocorrelations start to diverge after about 30 time units, 
and they also start to diverge from the autocorrelations of case (2) after about the same 
propagation time. Interestingly however, the least populated natural orbitals for logical 
coordinates Q^' 1,3 ' 3 ' 3 ' 3 and Q{ ' ' ' ' have populations of the order of 10~ 3 to 10~ 4 , and this 
is substantially improved to 10~ 4 to 10~ 5 for the non-homogeneous wavefunction. All other 
natural populations in the whole tree structure, for both wavefunctions, remain in the order 
of 10~ 4 or below. The spectra of cases (1) and (2) for the non homogeneous wavefunctions 
(Fig. [8H) agree to each other rather well. The differences in the autocorrelations are however 
reflected in a appearance of a few spurious lines with low intensity in the low energy part of 
the spectrum for case (1). 
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The analysis of case (2) indicates that when correlated DOF are close to each other in the 
tree structure, i.e. they are connected in the deeper layers of the tree, it is easy to achieve 
convergence by increasing the number of SPFs close to such DOF in the tree structure. 
Even most importantly, one can expect an early convergence to the correct trend, since even 
with not very good natural populations at low layers the results are already correct. Case 
(1), on the contrary, is an example of the wrong choice of tree structure. Even with natural 
populations smaller than case (2) for the homogeneous wavefunction, and populations of the 
order of 10 -4 in most layers, the results are different from case (2) and also differ between the 
two different wavefunctions tried for case (1). In case (1) a general and substantial increase 
of SPFs for the whole tree structure would probably be needed to obtain better converged 
results. The early convergence of the method with the number of SPFs appears to be 
damaged in case (1) in light of the difference in results of for the two different wavefunctions 
tried. 



B. Pyrazine 



We now turn our attention to the photophysics of pyrazine using the 24D second-order 
vibronic-coupling Hamiltonian of Raab et al. [6|. In pyrazine, the spectrum obtained from 
excitation to the second excited electronic state S2 presents a broad feature due to the fast 
decay of the system into the S\ electronic state through a conical intersection. Such a decay 
occurs during a few ten fs after photoexcitation. Although obtaining the right position 
and approximate shape of the broad band is relatively straightforward and several different 
approaches can claim to have achieved it, reproducing the fine details of the absorption 
spectrum is quite hard, since a good quality propagation up to relatively long times is needed. 
Therefore, the 24D model of pyrazine in Ref. [6] has often been used to benchmark different 



quantum dynamical approaches including semic 



method 



50( , the coupled co 



rerent state method 



Fourier-transform method 



assical methods 



51 



49[, the multiple-spawning 



52J, the matchin g-pu rsuit / split-operator 



53|], and the Gaussian-MCTDH method [54]. 



In the 24D simulations reported here, we use the same DVR and grid points as in the 



MCTDH simulations of Refs. jo] and 54|. These MCTDH calculations were based on 
mode combination scheme in which the 24 primitive coordinates were grouped into the 
8 combined modes Q x = [vi 0a ,z/ 6o ], Q 2 = [1/1, i/ 9o , z/ 8a ], Qa = [^2,m,^8b], Qa = [^4,^5,^3], 
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and Q 8 = [v 20b , u 16b , v u , u 7b ]. 
6[ and it is beyond the purpose 



Q5 — [^160, ^12, ^13], — [ v l9bi ^18&], Ql — [^180) ^14, v VAai v . 

A description of these vibrational modes is found in Ref. 
of this paper to reproduce their discussion here. All but one of the reported ML-MCTDH 
simulations are based on this set of combined modes, which are then grouped in the ML tree 
as depicted in Fig. |9j This makes the ML-MCTDH calculations numerically easier to com- 
pare to previous MCTDH calculations, for example when computing wavefunction overlaps, 
since the ML-MCTDH wavefunction can be easily expanded into the corresponding MCTDH 
form (if the expanded form is still of a manageable size) in a similar way as an MCTDH 
wavefunction can be expanded to a standard wavefunction on the primitive grid, Eq. ([I]). In 
the simulation named ML-MCTDH-1, however, we further divide modes Q7 and Qs, which 
contain coordinates with large primitive grids, in order to better exploit the capabilities of 
ML-MCTDH. We used Qi---Q 5 as above but: Q 6 = [u 19b ], Q 7 = [u 18b ], Q 8 = [v 18a ,v 14 \ : 
Q9 = [vwa^na], Q10 = Kefe], and Q n = [^206,^11,^76]- This division is accomplished by 
adding a further (a sixth) layer to the tree depicted in Fig. [9j 

FIGURE M AROUND HERE 

For non-adiabatic systems with more than one electronic state MCTDH is often used in its 



multi-set formulation 



3|-|5|, in which different sets of SPFs are used for each electronic DOF. 
All pyrazine MCTDH calculations discussed here were done using the multi-set formalism. 
In the multi-set formalism the SPFs of each state are optimal for that state, so that a smaller 
number of SPFs per state are needed than in a single-set calculation. The computational cost 
grows exponentially with the product of SPFs in each state, but linearly with the number of 
SPFs in different electronic states, which makes multi-set calculations often advantageous 
with respect to single-set ones for the same system. In ML-MCTDH, a multi-set formulation 
would be extremely cumbersome because one would end up with several tree structures 
to be specified, one for each electronic state, and a much more complex algorithm. On 
the other hand, the favorable scaling of ML-MCTDH with the number of DOF makes it 
unnecessary to use a multi-set ML formulation. In simulations of non-adiabatic systems 
with ML-MCTDH the electronic DOF is hence just another coordinate that indexes the 
electronic states, as in usual single-set MCTDH calculations, and the ML-MCTDH algorithm 
remains unchanged. Consequently, in constructing the ML tree, one is free to group the 
electronic DOF in any convenient way with other coordinates and set it in any level of the 
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tree where it may seem appropriate. For example, if only a subset of coordinates couple 
the different electronic states and the rest act as a bath, it may make sense to group the 
electronic DOF with the coordinates that control the non-adiabatic coupling, and then the 
system-bath separation can be placed one level above. In our simulations, the top layer 
contains coordinates Q V ib, q e i, where Q V ib groups the 24 vibrational modes of the system and 
q e i denotes the electronic DOF. Q vi f, is further divided into two groups of coordinates, one 
containing modes z/ 10a , z/ 6a , ui, u 9a , z/ 8a and the other one containing the rest of modes. Modes 
^lOa; v&a-, v \i ^9a define the "system" in 4D models of pyrazine |6|, while the rest have been 
usually termed "bath" modes. Therefore, the first branching of the tree after separating 
electronic and vibrational coordinates can be understood as a system-bath separation. The 
resulting system and bath coordinates are further divided until modes of a manageable size 
are reached, as seen in Fig. 

TABLE E AROUND HERE 

Table |T] presents some results for a set of ML-MCTDH simulations and three different 
MCTDH reference calculations. The calculation named here MCTDH-1 was reported in 



Ref. The MCTDH-2 result was an MCTDH reference calculation in Ref. [54j, while 
the MCTDH-3 calculation is a new reference result for this system and corresponds to the 
largest MCTDH calculation reported for pyrazine to date. All simulations were run up to a 
propagation time of 150 fs. 

FIGURE [TO] AROUND HERE 

FIGURE E] AROUND HERE 

Fig. HO] presents the absolute value of the autocorrelation function for the ML-MCTDH and 
MCTDH calculations in Table HI As in the HH cases above, the autocorrelation is computed 
using the relation for real initial states and symmetric Hamiltonians a(t) = (\I/*(t/2)|\I/(t/2)). 
In Fig. [TUb one sees that the ML-MCTDH-1 simulation presents severe deficiencies in |a(£)|, 
only following the general shape of the initial decay and the small revivals between 25 and 45 
fs. This simulation is however extremely fast, needing only 7 minutes to run and consuming 
an incredibly small amount of resources. Regarding peak positions and intensities, the 
general features of the spectrum in Fig. [TTk . although not too accurate, are there. The 
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details of the peaks below 2.1 eV are not reproduced and the high energy tail, above 2.5 eV, 
shows an artificial oscillatory behavior. But beside this, the main peak is amazingly well 
reproduced, considering the inexpensiveness of the calculation. The larger ML-MCTDH-2 
simulation, needing 33 hours (five times faster than the MCTDH-1 reference result), already 
presents the correct trends in the autocorrelation function up to 70 fs, resulting in a much 
improved spectrum (Figs. IT2"r and [T2"b). The autocorrelations of the two most accurate 
ML-MCTDH simulations present features that resemble closely those of the most accurate 
MCTDH result. In fact, when comparing the various structures of \a(t)\ in Fig. [TUb . the 
tendency of the MCTDH-1, -2, and -3 simulations is towards the ML-MCTDH-8 simulation, 
which seems to indicate that the best ML-MCTDH results are better converged than the 
best MCTDH run. This can be seen in the height of the structure between 30 and 40 fs, 
and also in the position of the three peak structures between 55 and 70 fs. When turning 
to the spectra in Fig. [T5J similar conclusions can be drawn, for example examining closely 
the peaks at 2.2 and 2.35 eV in Fig. [12b. It is remarkable that simulation ML-MCTDH-7, 
which seems to be already better converged than the best MCTDH result reported, required 
about 25% of runtime and 1% of memory (in terms of wavefunction size) than the latter. 

FIGURE [H AROUND HERE 

In terms of efficiency, the second column in Table [I] contains the wall-clock time of each 
calculation. All ML-MCTDH calculations were run on a single-processor, so CPU time 



equals in this case wall-clock time, and used the variable mean- field (VMF) j^] propagation 



scheme, 
code 



he MCTDH calculations were run using the shared-memory parallelized MCTDH 



55|,I5 

e. |a 



56|and 8 processors in parallel, using the constant mean-field (CMF) propagation 



scheme. [3j, |57j, |58| The speed-up factor of the MCTDH parallel computations is 2.9, 3.3, and 
3.7 for MCTDH-1, -2, and -3, respectively, and this factor has been already multiplied to the 
MCTDH wall-clock times, yielding the CPU time that would have been taken on a single 
CPU so that they can be readily compared to the ML-MCTDH values. These speed-up 
are reasonable but not particularly good, which is due to the relatively small Hamiltonian 
operator of the pyrazine model and better scalings have been obtained, e.g. for the Zundel 



cation H 5 2 



55 



56 



59J. All calculations were run on the same machine and CPU type 
(see caption of Table H]). All reported simulations, including the HH ones but except for 
ML-1, were run with a rather high integrator precision of 10 -7 . This was done to exclude 
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any numerical artifacts and to ensure that all discussed effects originate form various sizes 
of sets of SPFs. The fast but low-accuracy calculation ML-1 was run with an integrator 
precision of 10 -5 . Reducing the high integrator precision will speed-up the calculations by 
factors between 1.5 and 2 without introducing visible changes into the spectra. 

The third column in Table [I] contains the total number of time- dependent coefficients 
in the wavefunction representation. Two aspects are remarkable here. First, ML-MCTDH 
wavefunctions are much more compact than MCTDH ones. Second, the CPU time divided 
by the number of coefficients is about one order of magnitude larger in ML-MCTDH, this 
factor remaining quite stable along the series of calculations. Therefore, integrating each 
coefficient in ML-MCTDH is harder than in MCTDH due to the natural overhead of the 
ML-MCTDH algorithm. However, this is compensated by the much smaller number of 
coefficients in the wavefunction. Such scaling effects become more and more pronounced 
with the size of the system, as was seen for the HH case, and ML-MCTDH propagations 
seem to become more efficient than MCTDH ones in terms of CPU time for systems larger 
than about 20 DOFs. Simulations ML-MCTDH-2 to -8 could be made even more efficient 
further dividing the largest combined modes at the lowest layer in a similar way as it was 
done for ML-MCTDH- 1. Running on 8 processors, MCTDH-3 needed about 1000 hours to 
complete (unsealed wall-clock time), keeping up with the wall-clock time of the best ML- 
MCTDH calculations through parallelization. The MCTDH-3 calculation could in principle 
be made faster than ML-MCTDH using more hardware. However, for a molecular system 
doubling the size of pyrazine, a normal MCTDH calculation would be probably not doable, 
or at least not with a reasonable SPF basis size. On the contrary, that case would require an 
additional layer in a ML-MCTDH calculation, and the CPU cost would scale approximately 
linearly, rendering the simulation still doable. Regarding the different approaches tried on 
the pyrazine Hamiltonian over the years, ML-MCTDH offers by far the best quality/cost 
relation. The ML-MCTDH-2 simulation, which is already of a reasonable quality, takes 
about one day on a single CPU using little resources. The various ML-MCTDH simulations 
on the cheaper end yield spectra of a quality that other approaches hardly reach 
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or reach only by using a much larger amount of CPU and memory resources [53|, [54(. And 
again, accepting slightly larger deviations in the spectrum, an ML-MCTDH calculation on 
pyrazine takes only 7 minutes. 
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IV. SUMMARY AND CONCLUSION 



We presented the implementation of the ML-MCTDH approach into the Heidelberg 
MCTDH package, which is based on the recursive algorithm proposed by Manthe. A dis- 
cussion of the concept of multilayer was given, and the working equations were provided 
making emphasis on the key points of the approach. For a comprehensive derivation of the 



working scheme one should see Ref 



411 ]. The use of the separable terms in the Hamiltonian 



in the recursive construction of the /i-matrices and mean-fields at the different layers was 
discussed here with some detail. 

The implementation, numerical performance, and general features of the ML-MCTDH ap- 
proach have been tested by running simulations on the Henon-Heiles model and on pyrazine. 
For Henon-Heiles, simulations for the 6D, 18D and a large 1458D case have been reported. 
In the 6D and 18D cases, normal MCTDH calculations are performed for identical param- 
eter sets for comparison, and for the 6D case MCTDH is more efficient than ML-MCTDH. 
The situation already changes for the 18D simulations, in which ML-MCTDH outperforms 
MCTDH. With respect to the natural populations, which are often used as a convergence 
criterion in MCTDH calculations, it is observed that when increasing the number of SPFs at 
lower layers for a given number of SPFs at higher layers, the convergence of the higher layer 
becomes worse although the overall quality may improve. This is due to the fact that the 
dynamics at the higher layers becomes more complex after increasing the size of the basis 
below. Hence, with respect to convergence, all layers have to be monitored when changing 
the basis size in some other layer. 

For pyrazine, we test ML-MCTDH on the second-order vibronic-coupling Hamiltonian 
of Raab et. al, which has been used in testing several quantum dynamical methods. By 
properly choosing the tree scheme, we report a ML-MCTDH simulation that takes only 7 
minutes of CPU time in one processor and very few memory, and that recovers the main 
spectral features reasonably well. On the other hand, we report ML-MCTDH simulations 
which are converged to at least the same degree as the best reference MCTDH results 
available, taking only 25% of the CPU time and 1% of wavefunction storage space with 
respect to such reference results. As seen in the pyrazine case, ML-MCTDH allows for a 
consistent separation of the degrees of freedom considered as system and the ones considered 
as bath, providing the variationally optimal system-bath evolution for the particular selection 
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of the basis size and the layering scheme. By using ML-MCTDH for system-bath problems, 



schemes that simplify the system-bath ansatz 
subset of effective modes 



or schemes that reduce a larger bath to a 



61] . become unnecessary. 
ML-MCTDH is an efficient tool for model systems, which had already been demonstrated 
by works of Wang and Thoss. ML-MCTDH has been shown here to be a very powerful tool 
to treat more realistic molecular Hamiltonians as well. The fact that a tree structure needs to 
be chosen, often based on intuition or experience, and that convergence has to be monitored 
at the different layers, makes its use less straightforward than MCTDH. Therefore, some 
work in the direction of automating certain decisions regarding the tree structure will have 
to be carried on. Regarding efficiency, it will be interesting to develop good dedicated 
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58] for the ML case. And, of course, 



constant mean-field (CMF) propagation schemes 
parallelizing the ML-code will be an important step, as ML-MCTDH is for treating large 
systems. Beyond applications of ML-MCTDH to problems involving nuclear dynamics, for 
which it has certainly a large potential, new directions will have to be explored in the future. 
We believe that methods based on or taking advantage of the ML-MCTDH concept can be 
useful in a wide range of situations. For example, we envisage that in simulations of mixtures 
of different kinds of particles (different kinds of fermions and bosons, electrons and nuclei 
in molecules, etc.), the different types of particles can be separated and correlated to each 
other at the top layer, while the internal dynamics of each group takes place in the layers 
below. ML-MCTDH constitutes a consistent and powerful tool for the quantum-dynamical 
description of high dimensional molecular systems. This and other kinds of challenging 
applications remain the subject of future investigations. 
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TABLE I: Simulation parameters of the various MCTDH and ML-MCTDH 24D pyrazine calcu- 
lations. The second column contains the wall-clock time of each simulation. ML-MCTDH calcula- 
tions were run on a single CPU, so that the wall-clock time equals the CPU time. The MCTDH 
calculations were run on 8 CPUs using shared-memory parallelization 



55, 



561 ] . The speed-up factor 



of the 8-processor parallel pyrazine calculations is 2.9, 3.3 and 3.7 for the MCTDH-1, -2 and -3 
cases, respectively. The wall-clock times given for the MCTDH reference results are already scaled 
up by the corresponding speed-up factor and therefore reflect the time that such simulation would 
have taken on a single processor. Therefore they can be readily compared to the ML-MCTDH 
values. All simulations were run on the same machine and CPU type, namely Quad-Core AMD 
Opterons, processor type 2384 running at 2.7 GHz. The third column shows the total number of 
time-dependent coefficients propagated in each case. The fourth column contains, for ML-MCTDH 
simulations, the number of SPFs for each node of the tree according to the representation in Fig. [9j 
The parenthesis indicate that different N$ values were taken for each of the branches. The asterisk 
for the ML-1 case indicates that there is a further layer below the branches corresponding to 
(See text and Fig. [9] for details). For the MCTDH calculations, the fourth column contains in 
each parenthesis the number of SPFs for combined modes Qi to Q% for electronic states S± and 
S2, respectively. All reported simulations, except for ML-1, were run with a rather high integrator 
precision of 10~ 7 . (See text). ML-1 was run with 10~ 5 . Reducing the high integrator precision 
will speed-up the calculations by factors between 1.5 and 2 without introducing visiable changes 
into the spectra. 
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Figure Captions 

Figure [TJ Tree structures for the MCTDH and ML-MCTDH wavefunctions of the 6D 
Henon-Heiles simulations, (a) MCTDH wavefunction tree, in which the coordinates are 
combined in groups of two. A^ refers in this particular case to the number of primitive 
basis functions or grid points, (b) ML-MCTDH wavefunction tree. This tree is similar to 
the MCTDH tree, but the three combined modes have been separated by adding an extra 
layer. For the ML-MCTDH wavefunction A" 3 corresponds now the number of primitive 
basis functions and for A2 = A3 (and same Ai) the ML-MCTDH case becomes numerically 
identical to the MCTDH one. 

Figure [2} (color) (a) Absolute value of the autocorrelation function for various MCTDH 
and ML-MCTDH 6D Henon-Heiles simulations with a coupling parameter A = Ao- The 
numbers in parenthesis correspond to (Ni, N 2 ) in Fig. [TJ In the key, plain numbers indicate 
MCTDH calculations and the symbol ML designates ML-MCTDH calculations, (b) Differ- 
ence between the autocorrelations of the ML-MCTDH simulations and the (50,24) MCTDH 
calculation. 

Figure [3t (color) (a) Absolute value of the autocorrelation function for various MCTDH 
and ML-MCTDH 6D Henon-Heiles simulations with a coupling parameter A = 2Ao- The 
numbers in parenthesis correspond to (Ni,N 2 ) in Fig. [TJ (b) Difference between the auto- 
correlations of the ML-MCTDH simulations and the (50,24) MCTDH calculation. In the 
key, ML is for ML-MCTDH. 

Figure [H Tree structures for the MCTDH and ML-MCTDH wavefunctions of the 18D 
Henon-Heiles simulations. The MCTDH wavefunction (a) consists of six combined modes, 
each of them grouping three primitive coordinates. The ML-MCTDH wavefunction (b) 
starts dividing the system in three logical coordinates, and each of them is further divided 
in three combined modes. The resulting tree has three layers of TD coefficients. 

Figure [5j (color) (a) Absolute value of the autocorrelation function for various MCTDH 
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and ML-MCTDH 18D Henon-Heiles simulations with a coupling parameter A = Ao- For 
the MCTDH cases the number in parenthesis refers to Aq in Fig. H^. For the ML-MCTDH 
cases the number in parenthesis refers to Aq and in Fig. Ub, which are set equal in the 
reported simulations, (b) Detailed view of the last recurrence structure between 55 and 60 
time units. In the key, ML is for ML-MCTDH. 



Figure HO (color) (a) Absolute value of the autocorrelation function for various MCTDH 
and ML-MCTDH 18D Henon-Heiles simulations with a coupling parameter A = 2A . For 
the MCTDH cases the number in parenthesis refers to Aq in Fig. H^. For the ML-MCTDH 
cases the number in parenthesis refers to Aq and N2 in Fig. |4b, which are equal in the 
reported simulations, (b) Detailed view of the recurrence structures between 10 and 30 time 
units. In the key, ML is for ML-MCTDH. 



Figure [7J (color) (a) Spectra for various MCTDH and ML-MCTDH 18D Henon-Heiles 
simulations with a coupling parameter of A = 2\q. The numbers in parenthesis have the 
same meanings as in Fig. [6j (b) Detailed view of the three first peaks between 6.2 and 7.9 
energy units. In the key, ML is for ML-MCTDH. The spectra are obtained by a Fourier- 
transform of the autocorrelation function using a cos 2 filter {3, 5] to minimize spurious effects 
known as Gibbs phenomenon. 



Figure [8j (color) Absolute value of the autocorrelation of the ML-MCTDH 1458D Henon- 
Heiles simulations for time units 30 to 60. (a) Reference results without coupling term 
between the two displaced coordinates for case (1) (red dashed), in which the two initially 
displaced degrees of freedom belong to different logical coordinates at all layers, and case 
(2) (blue line), in which the two initially displaced degrees of freedom belong to the same 
logical coordinate at all layers, (b) For case (2), homogeneous wavefunction (red dashed) 
and non homogeneous wavefunction with more SPFs at low layers (blue line), (c) For case 
(1), homogeneous wavefunction (red dashed) and non homogeneous wavefunction (blue line). 
The pairs of autocorrelation functions presented in plots (b) and (c) are almost identical 
until about 15 time units, but the autocorrelation functions of case (1) and (2) start to differ 
already after 5 time units. These differences increase with time and can be inspected for 
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the time interval 30 to 60 by comparing figure (b) with (c). (d) Spectra of the calculations 
with the non homogeneous wavefunctions for cases (1) and (2). The spectra are obtained 
by a Fourier-transform of the autocorrelation function using a cos filter to minimize 

spurious effects known as Gibbs phenomenon. 

Figure [9t Tree structure used in most of the ML-MCTDH simulations of 24D pyrazine. 
The maximum depth of the tree is five layers, and the first one separates the 24 vibrational 
coordinates and the discrete electronic degree of freedom. The number of SPFs denoted N± 
need not be necessarily the same, and in fact in some of the ML-MCTDH simulations they 
are chosen to be different. The same is true for the three N&S, which can be different. For 
the fast ML-1 calculation a 6th layer was added, see text. 

Figure HOt (color) Absolute value of the autocorrelation function for the 24D pyrazine cal- 
culations. The numeration of the calculations is consistent with Table [B (a) Best MCTDH 
result and two best ML-MCTDH results from to 300 fs. (b) Detailed view of the recur- 
rences between 20 and 70 fs for all reported MCTDH simulations and the best ML-MCTDH 
simulation (b), and for the two worst and two best converged ML-MCTDH runs (c). In the 
key, ML is for ML-MCTDH. 

Figure [TTJ (color) (a) Spectra of the fastest ML-MCTDH 24D pyrazine calculation (2.2 x 
10 4 TD coeff., 7 min. of CPU) and the best reference MCTDH result (4.6 x 10 7 TD coeff., 
2901 hours of CPU), (b) Detailed view of the energy domain between 2.1 and 2.6 eV. The 
spectra are obtained by a Fourier-transform of the autocorrelation function using a cos 2 
filter 

aa 

to minimize spurious effects known as Gibbs phenomenon. The numeration of 
the calculations is consistent with Table [H In the key, ML is for ML-MCTDH. 

Figure 112b (color) (a) Spectra of the three reference MCTDH simulations and three of 
the ML-MCTDH simulations, among them the two best results, for 24D pyrazine. The 
spectra are obtained by a Fourier-transform of the autocorrelation function using a cos 2 
filter 

aa 

to minimize spurious effects known as Gibbs phenomenon. The numeration of 
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the calculations is consistent with Table HI For the energy domain between 2.1 and 2.6 eV: 

(b) The three reference MCTDH calculations compared to the best ML-MCTDH result and 

(c) ML-MCTDH simulations 2, 7 and 8. In the key, ML is for ML-MCTDH. 
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FIG. 2: Vendrell and Meyer, Journal of Chemical Physics 
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FIG. 3: Vendrell and Meyer, Journal of Chemical Physics 
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FIG. 4: Vendrell and Meyer, Journal of Chemical Physics 
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FIG. 5: Vendrell and Meyer, Journal of Chemical Physics 
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FIG. 7: Vendrell and Meyer, Journal of Chemical Physics 
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